Answer the following questions:
In this activity, you will:
O’Sullivan D and Unwin D (2010) Geographic Information Analysis, 2nd Edition, Chapter 7. John Wiley & Sons: New Jersey.
For this activity you will need the following:
Hamilton CMA CT.It is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:
rm(list = ls())
Note that ls() lists all objects currently on the worspace.
Load the libraries you will use in this activity. In addition to tidyverse, you will need spdep, a package designed for the analysis of spatial data (you can learn about spdep here and here):
library(tidyverse)
-- Attaching packages --------------------------------- tidyverse 1.2.1 --
v ggplot2 2.2.1 v purrr 0.2.4
v tibble 1.4.2 v dplyr 0.7.4
v tidyr 0.8.0 v stringr 1.2.0
v readr 1.1.1 v forcats 0.2.0
-- Conflicts ------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(rgdal)
Loading required package: sp
rgdal: version: 1.2-16, (SVN revision 701)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
Path to GDAL shared files: C:/Users/Antonio/Documents/R/win-library/3.4/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: C:/Users/Antonio/Documents/R/win-library/3.4/rgdal/proj
Linking to sp version: 1.2-7
library(broom)
library(spdep)
Loading required package: Matrix
Attaching package: <U+393C><U+3E31>Matrix<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:
expand
Loading required package: spData
Begin by loading the data that you will use in this activity:
Hamilton_CT <- readOGR(".", layer = "Hamilton CMA CT", integer64 = "allow.loss")
OGR data source with driver: ESRI Shapefile
Source: ".", layer: "Hamilton CMA CT"
with 188 features
It has 255 fields
Integer64 fields read as signed 32-bit integers: ID POPULATION PRIVATE_DW OCCUPIED_D ALL_AGES AGE_4 AGE_5_TO_9 AGE_10_TO_ AGE_15_TO_ AGE_15 AGE_16 AGE_17 AGE_18 AGE_19 AGE_20_TO_ AGE_25_TO_ AGE_30_TO_ AGE_35_TO_ AGE_40_TO_ AGE_45_TO_ AGE_50_TO_ AGE_55_TO_ AGE_60_TO_ AGE_65_TO_ AGE_70_TO_ AGE_75_TO_ AGE_80_TO_ AGE_85 MEDIAN_AGE MALE_ALL_A MALE_4 MALE_5_TO_ MALE_10_TO MALE_15_TO MALE_15 MALE_16 MALE_17 MALE_18 MALE_19 MALE_20_TO MALE_25_TO MALE_30_TO MALE_35_TO MALE_40_TO MALE_45_TO MALE_50_TO MALE_55_TO MALE_60_TO MALE_65_TO MALE_70_TO MALE_75_TO MALE_80_TO MALE_85 MALE_MEDIA FEMALE_ALL FEMALE_4 FEMALE_5_T FEMALE_10_ FEMALE_15_ FEMALE_15 FEMALE_16 FEMALE_17 FEMALE_18 FEMALE_19 FEMALE_20_ FEMALE_25_ FEMALE_30_ FEMALE_35_ FEMALE_40_ FEMALE_45_ FEMALE_50_ FEMALE_55_ FEMALE_60_ FEMALE_65_ FEMALE_70_ FEMALE_75_ FEMALE_80_ FEMALE_85 FEMALE_MED MARRIED_AG MARRIED_OR MARRIED COMMON_LAW UNMARRIED SINGLE SEPARATED DIVORCED WIDOWED MARRIED_A1 MARRIED_O1 MARRIED_M COMMON_LA1 UNMARRIED_ SINGLE_M SEPARATED_ DIVORCED_M WIDOWED_M MARRIED_A2 MARRIED_O2 MARRIED_F COMMON_LA2 UNMARRIED1 SINGLE_F SEPARATED1 DIVORCED_F WIDOWED_F FAMILIES_I FAMILY_SIZ FAMILY_SI1 FAMILY_SI2 FAMILY_SI3 COUPLE_FAM COUPLE_MAR COUPLE_MA1 COUPLE_MA2 COUPLE_MA3 COUPLE_MA4 COUPLE_MA5 COUPLE_COM COUPLE_CO1 COUPLE_CO2 COUPLE_CO3 COUPLE_CO4 COUPLE_CO5 SINGLE_PAR SINGLE_PA1 SINGLE_PA2 SINGLE_PA3 SINGLE_PA4 SINGLE_PA5 SINGLE_PA6 SINGLE_PA7 SINGLE_PA8 CHILDREN_F CHILDREN_1 CHILDREN_2 CHILDREN_3 CHILDREN_4 CHILDREN_5 POPULATIO1 POPULATIO2 POPULATIO3 POPULATIO4 POPULATIO5 POPULATIO6 POPULATIO7 POPULATIO8 POPULATIO9 POPULATI10 POPULATI11 POPULATI12 PRIVATE_HO PRIVATE_HH PRIVATE_H1 PRIVATE_H2 PRIVATE_H3 PRIVATE_H4 PRIVATE_H5 PRIVATE_H6 PRIVATE_H7 PRIVATE_H8 PRIVATE_H9 PRIVATE_10 PRIVATE_11 PRIVATE_12 PRIVATE_13 PRIVATE_14 PRIVATE_15 OCC_PRIVAT OCC_PRIVA1 OCC_PRIVA2 OCC_PRIVA3 OCC_PRIVA4 OCC_PRIVA5 OCC_PRIVA6 OCC_PRIVA7 OCC_PRIVA8 OCC_PRIVA9 PRIVATE_16 PRIVATE_17 PRIVATE_18 PRIVATE_19 PRIVATE_20 PRIVATE_21 PRIVATE_22 PRIVATE_23 NATIVE_LAN NATIVE_LA1 NATIVE_LA2 NATIVE_LA3 NATIVE_LA4 NATIVE_LA5 NATIVE_LA6 NATIVE_LA7 NATIVE_LA8 NATIVE_LA9 NATIVE_L10 NATIVE_L11 NATIVE_L12 NATIVE_L13 NATIVE_L14 NATIVE_L15 NATIVE_L16 NATIVE_L17 NATIVE_L18 NATIVE_L19 NATIVE_L20 NATIVE_L21 NATIVE_L22 NATIVE_L23 NATIVE_L24 NATIVE_L25 NATIVE_L26 NATIVE_L27 NATIVE_L28 NATIVE_L29 NATIVE_L30 NATIVE_L31 NATIVE_L32 NATIVE_L33 NATIVE_L34 NATIVE_L35 NATIVE_L36 NATIVE_L37 NATIVE_L38 NATIVE_L39 NATIVE_L40 NATIVE_L41 NATIVE_L42 NATIVE_L43 NATIVE_L44 NATIVE_L45 NATIVE_L46 NATIVE_L47 NATIVE_L48 NATIVE_L49 NATIVE_L50 NATIVE_L51 NATIVE_L52 NATIVE_L53 NATIVE_L54 NATIVE_L55 NATIVE_L56 NATIVE_L57 NATIVE_L58 NATIVE_L59
You can obtain new (calculated) variables as follows. For instance, to obtain the proportion of residents who are between 20 and 34 years old, and between 35 and 49:
Hamilton_CT@data <- mutate(Hamilton_CT@data, Prop20to34 = (AGE_20_TO_ + AGE_25_TO_ + AGE_30_TO_)/POPULATION, Prop35to49 = (AGE_35_TO_ + AGE_40_TO_ + AGE_45_TO_)/POPULATION)
Example: Proportion of residents who are between 20 and 34 years old, and between 35 and 49:
Hamilton_CT@data <- dplyr::transmute(Hamilton_CT@data, AREA = AREA, TRACT = TRACT,
POPULATION = POPULATION,
POP20to34 = (AGE_20_TO_ + AGE_25_TO_ + AGE_30_TO_),
Prop20to34 = POP20to34/POPULATION,
POP35to49 = (AGE_35_TO_ + AGE_40_TO_ + AGE_45_TO_),
Prop35to49 = POP35to49/POPULATION,
POP50to64 = (AGE_50_TO_ + AGE_55_TO_ + AGE_60_TO_),
Prop50to64 = POP50to64/POPULATION,
POP65Plus = (AGE_65_TO_ + AGE_70_TO_ + AGE_75_TO_ + AGE_80_TO_ + AGE_85),
Prop65Plus = POP65Plus/POPULATION)
This is a SpatialPolygonDataFrame. Convert to a dataframe (“tidy” it) for plotting using ggplot2:
Hamilton_CT.t <- tidy(Hamilton_CT, region = "TRACT")
Hamilton_CT.t <- rename(Hamilton_CT.t, TRACT = id)
Rejoin the data:
Hamilton_CT.t <- left_join(Hamilton_CT.t, Hamilton_CT@data, by = "TRACT")
Column `TRACT` joining character vector and factor, coercing into character vector
This is the function to create local Moran maps:
localmoran.map <- function(spat_pol = spat_pol, listw = listw, VAR = VAR, ID = ID){
require(tidyverse)
require(broom)
require(spdep)
require(plotly)
spat_pol@data <- data.frame(ID = ID, VAR = VAR)
spat_pol.t <- broom::tidy(spat_pol, region = "ID")
spat_pol.t <- dplyr::rename(spat_pol.t, ID = id)
spat_pol.t <- dplyr::left_join(spat_pol.t, spat_pol@data, by = "ID")
df_msc <- transmute(spat_pol@data,
ID = ID,
Z = (VAR-mean(VAR)) / var(VAR),
SMA = lag.listw(listw, Z),
Type = factor(ifelse(Z < 0 & SMA < 0, "LL",
ifelse(Z > 0 & SMA > 0, "HH", "HL/LH"))))
local_I <- localmoran(spat_pol$VAR, listw)
spat_pol.t <- left_join(spat_pol.t,
data.frame(ID = spat_pol$ID, local_I))
spat_pol.t <- rename(spat_pol.t, p.val = Pr.z...0.)
spat_pol.t <- left_join(spat_pol.t,
df_msc)
map <- ggplot(data = spat_pol.t,
aes(x = long, y = lat, group = group,
p.val = p.val, VAR = VAR)) +
geom_polygon(aes(fill = Type, color = p.val < 0.05)) +
scale_fill_brewer(palette = "RdBu") +
scale_color_manual(values = c(NA, "Black") ) +
labs(color = "Prob < 0.05") +
coord_equal() +
theme(legend.title = element_blank())
ggplotly(map, tooltip = c("p.val", "VAR"))
}
This is function is used to create \(G_i^*\) maps:
gistar.map <- function(spat_pol = spat_pol, listw = listw, VAR = VAR, ID = ID){
require(tidyverse)
require(broom)
require(spdep)
require(plotly)
spat_pol@data <- data.frame(ID = ID, VAR = VAR)
spat_pol.t <- broom::tidy(spat_pol, region = "ID")
spat_pol.t <- dplyr::rename(spat_pol.t, ID = id)
spat_pol.t <- dplyr::left_join(spat_pol.t, spat_pol@data, by = "ID")
df.lg <- localG(VAR, listw)
df.lg <- as.numeric(df.lg)
df.lg <- data.frame(Gstar = df.lg, p.val = 2 * pnorm(abs(df.lg), lower.tail = FALSE))
df.lg <- mutate(df.lg,
Type = factor(ifelse(Gstar < 0 & p.val <= 0.05, "Low Concentration",
ifelse(Gstar > 0 & p.val <= 0.05, "High Concentration", "Not Signicant"))))
spat_pol.t <- left_join(spat_pol.t,
data.frame(ID = spat_pol$ID, df.lg))
map <- ggplot(data = spat_pol.t,
aes(x = long, y = lat, group = group,
p.val = p.val, VAR = VAR)) +
geom_polygon(aes(fill = Type, color = p.val < 0.05)) +
scale_fill_brewer(palette = "RdBu") +
scale_color_manual(values = c(NA, "Black") ) +
labs(color = "Prob < 0.05") +
coord_equal() +
theme(legend.title = element_blank())
ggplotly(map, tooltip = c("p.val", "VAR"))
}
Create spatial weights.
Hamilton_CT.w <- nb2listw(poly2nb(pl = Hamilton_CT))
Hamilton_CT.3knb <- Hamilton_CT %>% coordinates() %>% dnearneigh(d1 = 0, d2 = 3, longlat = TRUE)
Hamilton_CT.3kw <- nb2listw(include.self(Hamilton_CT.3knb), style = "B")
localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)
Loading required package: plotly
Attaching package: <U+393C><U+3E31>plotly<U+393C><U+3E32>
The following object is masked from <U+393C><U+3E31>package:ggplot2<U+393C><U+3E32>:
last_plot
The following object is masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:
filter
The following object is masked from <U+393C><U+3E31>package:graphics<U+393C><U+3E32>:
layout
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$Prop20to34, Hamilton_CT$TRACT)
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
gistar.map(Hamilton_CT, Hamilton_CT.3kw, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
gistar.map(Hamilton_CT, Hamilton_CT.3kw, Hamilton_CT$Prop20to34, Hamilton_CT$TRACT)
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34/Hamilton_CT$AREA, Hamilton_CT$TRACT)
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`